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Abstract In this paper we revisit and adapt to viral evolution an approach based on the the- 
ory of branching process advanced by Demetrius, Schuster and Sigmund ("Polynucleotide 
evolution and branching processes". Bull. Math. Biol. 46 (1985) 239-262), in their study of 
polynucleotide evolution. By taking into account beneficial effects we obtain a non-trivial 
multivariate generalization of their single-type branching process model. Perturbative tech- 
niques allows us to obtain analytical asymptotic expressions for the main global parameters 
of the model which lead to the following rigorous results: (i) a new criterion for "no sure 
extinction", (ii) a generalization and proof, for this particular class of models, of the lethal 
mutagenesis criterion proposed by Bull, Sanjuan and Wilke ("Theory of lethal mutagenesis 
Q ■ for viruses", J. Virology 18 (2007) 2930-2939), (iii) a new proposal for the notion of relax- 

ation time with a quantitative prescription for its evaluation, (iv) the quantitative description 
of the evolution of the expected values in in four distinct "stages": extinction threshold, 
O^' lethal mutagenesis, stationary "equilibrium" and transient. Finally, based on these quantita- 

tive results we are able to draw some qualitative conclusions. 
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1 Introduction 

RNA viruses exhibit a pronounced genetic diversity dPomingo and et'al. This vari- 

abiUty allows RNA vi ms to better adapt to en vironmental challenges as represented by host 
and therapy pressures dPomingo et all Il998h . Due to the lack of a proofreading activity 
of viral RNA polyme rases (average e r ror in corporation rate of order lO^'* per nucleotide, 
per replication cycle dSteinhauer et all Il992h ). short generation times and huge population 
numbers, RNA viral populations may be viewed as a collection of particles bearing mutant 
genomes. 

As a consequence of high mutation rates, frequencies of mutants depend not only on 
their level of adaptation but on the proba bility of being faithfully replicated during viral 
genomic RNA synthesis. Several studies dPomingo et all 1 19981 : lElena and Saniuail l2005t) 
have looked at viral diversification processes as a contributing cause of disease progression 
and of therapeutic strategies shortcomings including vaccine trials. It has become important 
to understand the process by which virus acquire diversity and the dynamics and fluctuations 
of this diversity in time. However, understanding viral evolution in vivo has proven to be a 
very unwieldy accomplishment due to the huge number of variables involved in the interplay 
between viruses and their hosts. 

For the last 30 years, quasi-species theory has provided a population-based framework 
to model RNA viral evolution. Quasi-species theory is a mathema tical framewo rk that was 
initially formulated to explain the evolution of macromolecules dEigenI , Il97ll) . More re- 
cently, it has been used to describe the evolutionary dynamics of RNA viruses. However, 
quasi-species theory is based on deterministic equations and so has a number of serious 
drawbacks when applied to realistic experimental systems since there are important sources 
of stochasticity that must be taken into account. 

Traditionally, in an effort to make the viral evolution process more palpable, several 
groups have addressed this subject from different points of view. There is a substantial 
amount of publications that studied virus populations during their evolution in experimen- 
tal settings, for instance, cell cultures (|Escaimfs et al, 2006), by challenging the virus with 



population bottlene cks dOiosnegros et all 2010bl al). or the introduction of antiviral drags 
dCrottv et a illooil), including mutagens, or another competing viral population. Several 
other groups have studied the proce ss of viral evolution away from the bench, using com- 
putational and mathematical tools dWilke et all 1200 ll ; iBull et all |2007[) . The main lesson 
learned from these efforts is that in order to escape from oversimplifying the interplay be- 
tween virus and hosts, one needs to take into account a few hard rales based on the ex- 
perimental data that has been generated by the community of investigators addressing viral 
evolution. 

In this paper we revisit and a dapt to viral evo l ution an approach based on the theory 
of branching process advanced by jPemetrius et a il dl985l) in their study of polynucleotide 
evolution. In fact, almost any stochastic model of asexual replication has an underlying 
branching process, which remains implicit and undefined in most of the studies in this sub- 
ject. For ins tance, (Manrabia et al', 2003:'Aguirre a nd ManrabiaL lioosi : I Aguirre et alll2009l : 
ICuesta etaii , ,201 UCap itan et al., .201 1 ; Cuesta, 201 ih have investigated probabilistic models 
in the form of a linear mean field approximation without explicitly defining the underly- 
ing stochastic process modeling the microscopic dynamics of particle replication. However, 
there are a few of them (Demetrius et al, 1985 ; Hermisson et al, 2002; Wilke, 2003) that take 
a different path and explicitly define the stochastic process in order to bring the mathemat- 
ical theory of branching processes to bear. This attitude has some virtues, since it provides 
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powerful tools, that have been perfected in the past several decades, allowing one to extract 
quantitative results on a rigorous basis in clear conc eptual framewor k. 

In the present w ork we pursue the same path as dOemetrius et all . [l985l : lHermisson et all 
l2002l : IWilkel[2003h , by putting more emphasis on the study of an explicit generating function 
of the underlying branching process and applying the mathematical tools as a means to get 
new insights and perspectives on the dynamics of evolving vims populations. Before setting 
up our mathematical framework let us give a broad overview of the scenery that we concieve 
as the context of the models described here. 

Viral infections may result in different clinical outcomes which are consequences of 
interactions between the virus and its host. The constant arms race between hosts and virus 
can either eliminate or perpetuate viral infections. Some viral diseases are self limited ending 
with the eradication of the infecting virus while others lead towards persistent infections. 

During its adaptation to new hosts virus populations will show an intense fluctuation in 
size and structure in response to the selective pressures imposed by the host defenses. Soon 
after the infection vims particles start to colonize the new environment and will seek within 
the host for the cell types that will better support the production of an exuberant progeny. 
At this moment and with the total absence of an immune response the viral population may 
experience an exponential growth. This growth is clinically reflected by the appearance of 
acute symptoms and the emergence of a peak of viremia detected in the peripheral blood. 
The emergence of selective pressures represented by the combination of components present 
in the innate and adaptive responses of the immune system will cause a drop in the viral load 
to levels well below the initial viremia value. At this point the vims and host defenses may 
reach a stationary "equilibrium " referred to as the viral set point in some chronic infections 
as the one caused by the human immunodeficiency vims (HIV). On the other hand, if the 
burden imposed by the host pressures overcomes the vims survival capacity, the viral popu- 
lation will become extinct. The initial exponential growth and subsequent population size re- 
duction experienced by the incoming vims can be referred to as the transient phase of a viral 
growth curve. The transient phase will represent the vims population fluctuations during the 
period from its introduction in the new host until the moment the viral set point is reached. 
In persistent infections without major oscillations, in the host selective pressures or in vims 
adaptability, vims populations will show an asymptotic behavior and will perpetuate in a 
stationary equilibrium regime. In vims with genome molecules represented by RNA which 
are able to tolerate high mutational rates, the stationary equilibrium is actually a balance 
between the two major forces in vims evolution: mutation and selection. As a matter of fact, 
the transient phase of a vims growth has been also called the recovery time as an allusion 
to vims populations recovering from bottleneck passages that correspond to drastic popu- 
lation size reductions. When selective pressures against the incoming vims are too strong 
for the vims to couple with them the stationary phase of the vims growth is never attained 
and the asymptotic behavior is absent. However, changes that will cause a recmdescence in 
the host selective pressures may dislocate vims from the stationary equilibrium towards its 
extinction. Therapeutic strategies as antiviral dmgs and therapeutic vaccines may precipi- 
tate vims populations in regions of instability represented by hostile adaptations marked by 
subsequent drops in vims replication capacity and progeny sizes. This phase of viral growth 
named here as extinction threshold represents a moment of uncertainty when the vims ex- 
tinction is strongly possible but unpredictable. Moving forward, when vims populations are 
unable to counteract the deleterious effects of all natural and chemical pressures combined, 
and average vims progeny drops to a value bellow one viable particle, the process called the 
lethal mutagenesis is started. At this moment the vims will become extinct almost certainly. 
Vims extinctions can be seen during self limited infections even before an asymptotic behav- 
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ior is reached or when chronic infections are cleared due to optimization of host pressures 
or to the addition of antiviral chemicals. 

Once the mathematical framework is established, the scenary described above raises 
several questions about the evolution of a viral population within its host can be posed and 
answered in a precise way. In this work we shall address the following questions: 

(i) What is the average number of particles in the population in each generation? 

(ii) What is the probability of extinction of the population within its host? 

(iii) How does the "transition" to extinction happens? 

(iv) When the population does not become extinct, what is the asymptotic behaviour? 

(v) How long it takes for the population to reach its asymptotic state? 

In order to answer these questions we propose a non-trivi al multivariate generalization of 
the single-type branching process of iDemetrius et a l([l98l. 

By employing perturbative techniques we are able to obtain analytical expressions for 
the main global parameter of the model, namely the average growth rate, represented here by 
th e malthusian par ameter of the branching process. Our results corroborate with the study 
of bull eTa 1 (l2007l^ by showing that the sufficient condition for lethal mutagenesis involves 
mutational and ecological aspects. They arrived at a conjectural criteria for lethal mutage- 
nesis by a heuristic and intuitive approach of great general applicability. For the particular 
class of models considered here it is possible to state and rigorously prove an extinction cri- 
te rion which is, un der a proper interpertation, equivalente to the lethal mutagenesis criterion 
of lBull eta il( l2007h . In other words the branching process formalism can provide a new and 
interesting perspective on this problem. We obtain a new criterion for non-extinction of a 
viral population and describe four distinct "stages" of evolution of a RNA vims populations: 
extinction threshold, lethal mutagenesis, stationary "equilibrium" and decay of the temporal 
auto-correlation. 



2 Phenotypic Models for Viral Evolution 

In this section we describe a model for viral evolution that is naturally represented by a 
multivariate branching stoch astic process generalizi ng, in a non-trivial way, the single-type 
branching process studied bv lPemetrius et al|( 1985h. For t he sak e of motivation we start by 



recalling a probabilistic model introduced bv lLazaroet"ail ( l2002l) . 

We interpret the notion of mutation probability as a general effect of probabilistic na- 
ture with direct impact on the replication capability of individual viral particles, considered 
here as a measure of the particle's fitness characterizing its phenotype. This effect is sum- 
marized by the definition of a stat ionary probability d i stribu tion which is used to set up a 
Galton- Watson branching process d Watson and GaltonL [l874l) for the temporal evolution of 
the viral population. This probability distribution gives appropriate parameters to classify 
the asymptotic behavior of the viral population and to describe some of the non-equilibrium 
properties of the model. 

In other publications on the same subject the concept of mutation is extensively used 
as the cause of replication capacity change. Understanding that those changes constitute 
an observable output due to many different factors (of genetic and non-genetic nature), we 
prefer to use the general term "effect" over the replication capacity to characterize the three 
possible changes (deleterious, beneficial and neutral) that may happen with the viral particle 
when it replicates. 
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2.1 The Probabilistic Model 

A number of viral infections starts with the transmission of a relatively small number of 
viral particles from one organism to another one. The initial viral population starts replicat- 
ing constrained by the unavoidable interaction with the host organism and evolves in time 
towards an eventual equilibrium. Each particle composing the population replicates in the 
cellular context that may differ from cell to cell. Moreover each particle has different replica- 
tion capabilities due to the natural genomic diversity found in viral populations in general. 
Therefore, it is reasonable to consider the viral population as a set of particles divided in 
groups of different replication capabilities measured in terms of the number of particles that 
one particle can produce. The models considered here do not take into account any informa- 
tion about the genomic diversity of any replicating class and therefore it should be classified 
as phenotypic model. 

Consider that the whole set of particles composing the viral population replicates at the 
same time in such a way that the evolution of the population is described as a succession 
of discrete viral generations. This assumption crucially depends on the clear definition of 
the time needed for a particle to replicate, referred by virologists as generation time. As 
it depends on the cell environment it is clear that this time period may vary from particle 
to particle, replicating in different cells, in such a way that a meaningful concept is a dis- 
tribution of replication times with a possible well defined mean value. The dispersion of 
the replication times can be considered small if we restrict ourselves to homogeneous cell 
populations. Under these conditions, one may consider that no particle can be part of two 
successive generations, that is, the generations are discrete and non-overlapping. 

Suppose that a population of viruses that start evolving from an initial set of particles, 
which is partitioned into classes according to their mean replication capacity, that is, for 
each integer r = 0, . . . ,R there is a class of particles labeled by r and a random variable 
assuming non-negative integer values whose probability distribution ti-{k) satisfies E(fr) = r 
and?o(^) = SkQ. 

Inasmuch as the process of replication is controlled by chemical reactions involving 
specific enzymes and the template, it is reasonable to assume a mean bounded replication 
capacity per particle that is possibly typical for each specific virus. Hence there is a maxi- 
mum mean replication capacity R imposed by the natural limiting conditions under which 
any particle of the population replicates. 

In the process of replication of a viral particle errors may occur at each replication 
cycle in the form of point mutations with possible impact on the replication capacity of 
the progeny particles. Due to the intrinsic stochastic component of chemical reactions it is 
natural to treat this point mutational cause as probabilistic. Another possible cause of change 
in the replication capability in the viral offspring is clearly related to the cellular environment 
where the replication process takes place. As a result the time evolution of viral populations 
should be viewed as a physical process strongly influenced by stochasticity. Therefore, the 
combined action of genetic and non-genetic causes may produce basically three types of 
replicative effects with associated probabilities, at the particles scale, applicable to every 
single replication event: 

- deleterious effect d: the mean replication capacity of the copied particle decreases by 
one. Note that when the particle has capacity of replication equal to it will not produce 
any copy of itself on the average. 
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- beneficial effect b: the replication capacity of the copy increases by one. If the mean 
replication capacity is already the maximum allowed then the mean replication capacity 
of the copies will stay the same. 

- neutral effect c: the mean replication capacity of the copies remain the same as the mean 
replication capacity of the parental particle. 

The only constraint these numbers should satisfy is b + c + d = l.ln the case of in vitro ex- 
periments with homogeneous cell populations the parameters c, d and b may be considered 
essentially as mutation probabilities. 

The assumption that b is very small when compared with d and c is justified by sev- 
eral experimental results. The frequencies between beneficial, deleterious and neutral mu- 



iMiralles et all 19991: Imhof and Schlotterei 


{ 200 ll Keiehtlev and LvnchL 2003 


; OrJ, 2003 


Saniuan et all 20041; Carrasco et all 12007; 


Evre- Walker and Keiahtlevl |2007|; 


Parera et al 


2007; Rokvtaetall 


20081). Taking their results together, it is reasonable to conclude that 



beneficial mutations could be as low as 1000 less frequent than either neutral or deleteri- 
ous mutations. As a result, the viral population would be submitted to a large number of 
successive deleterious and neutral changes and a comparatively small number of beneficial 
changes. The particular case where there are no beneficial effects in time (b = 0) is interest- 
ing not only because of biological reasons but also due to a mathematical property, namely, 
the spectral problem associated to its mean matrix is "completely solvable". This property 
will be crucial for us to implement the perturbative expansions in the general case where 
fo^O. 



2.2 The Multitype Branching Process 

Instead of working directly with a probabilistic model we will introduce a new generat- 
ing function, which is a non-trivial multivari ate generalization of the single-type branching 
process proposed by Demetrius et al ('1985'). For full presentations of the theory of branch- 
ing processes see (Harris, 1963; AthreyAand Ney, 1972; Kimrnel and Axelrod, 200i). Brief 
accounts of the theory branching p rocess, high hghting the main resu lts we need here, are 
provided by dOemetrius et alLll985l Sec. 2) and I Antoneli et all (l2013ah . 

In order to find our probability generating function we observe that the replication pro- 
cess is simply a Bernoulli trial with three possible outcomes: a newly produced particle may 
have endured a deleterious, neutral or beneficial eff ect. Therefore, the offspring distribution 
should be a trinomial distribution (see lFelleil(ll968h ) and the probability generating function 
of the phenotypic model is / = (/o,/i , . . . with components 

Mzq,Zu...,Zr) = ts{k){dz.s-l+CZs + bZs+l)''- (1) 

with "consistency conditions" /o(zo,Zi, ••■ = 1 andz«+i = zr- 

Note it is trivial to further generalize the model in order to include "higher order replica- 
tive effects" that changes the mean replication capacity by ±2, ±3, . . . and thus replacing 
the trinomial distribution by a multinomial distribution. It is also worth to mention that there 
are other possible variations of these models that share the same essential properties and are 
more adequate in different contexts: 
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With Zero Class: In this version the model naturaly fits the symmetry of the binomial dis- 
tribution since it contains the particles of class r = which are generated by the particles 
from class r = I. 

Without Zero Class: In this variation, the particle class is omitted and thus the probabil- 
ity generating function has R variables and R components: set the variable zo = 1 and 
omit the first component /q. Particles of class r = 1 undergoing a deleterious change are 
eliminated in the next generation. In this formulation the model is positively regular. 

With the convention described above, it is easy to see that in the one-dimensional case 
with b = one obtains the following generating function 

fiz) = £rW((l-c)+cz)^ = f^t(k){l-c{l-z))'. (2) 

k=0 A=0 

This i s exactly the generating function of the single-type model proposed by jPemetrius et all 
1 19851 . p. 255, eq. (49)) for the evolution of polynucleotides. In their formulation, c = is 
the probability that a given copy of a polynucleotide is exact, where the polymer has chain 
length of V nucleotides and there is a fixed probability p of copying a single nucleotide cor- 
rectly. The replication distribution t provides the number of copies a polynucleotide yields 
before it is degraded by hydrolyses. 



2.3 Evolution of the Mean Values 



We shall introduce the notation Z; 



,Q — 1 for the condition Zq = e,., which is the initial 
population consisting of one particle of class r and no particles of other classes. Thus 
P(Zi = i\ZQ = 1) = A basic assumption in the theory of branching processes is that 
all the first moments are finite and that they are not all zero. Then one may consider the 
mean evolution matrix or the matrix of first moments M = {Mij} which describes how the 
averages (Z„) of the sub-populations of particles in each replication class evolves in time: 



{Z„)=M"Zo or (Z„)=M(Z„_i). 
The mean matrix M = {Mij} is defined as 



(3) 



M,v = E(Zi|Z^ = l) 



or 



«„^|i(i), 



fori,; = 0,...,7?andl = (1,1,...,!). 

In order to compute the mean matrix of ([!} we observe that only the partial derivative 
of a component fr with respect to Zr-i, Zr and Zr+i are non-trivial, the others are zero. 
Therefore, the mean matrix is 



M - 



/Od 
c2d 
Ob 2c 3d 
2b 3c 

VOO 








{R-l)c 
(R-\)bR{ 








Rd 

c + b)J 



(4) 



Thus showing that our generating function indeed provides a non-trivial (i.e., non-trianglar) 
generalization of (|2}. Interestingly, the mean matrix M provided by this class of models is 
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a tridiagonal matrix, which is an ubi quitous type of matrix appearing in several fields rang- 
ing fro m statistical signal processing Gray (l2006h, infor mation theory iGrenander and Szegg 
( Il958l ). lattice dynamical systems dAntoneli et all l2005l) . 

The mean matrix obtained ab oye coincide wit h some o f the matrices in the line a r mea n 
field approx i matio n studied in Manmbia et all l2003h: Aguirre and Marrrubial ( 1200 8l) , 
lAguirre et all (l2009l) ; ICuesta et all (|201 lh : ICapitan et all(l201 ihJCuestd ( |201 ll) . Howeyer. it is 
important to stress that these works consider that the total population Z„ eyolyes by multi- 
plication by the mean matrix, while here it is the evolution of the average value (Z„) that 
is described by the mean matrix. However, the mean matrix M depends on the probability 
distributions tr{k) only through its mean value r. Therefore, there are infinitely many distinct 
branching processes with the same mean matrix and the only way to tell them apart is by 
looking at their fluctuations or second moment properties. 

Let p {M) denote the spectral radius of the mean matrix M, that is, if Ai ^ • • • Xr are 
the eigenvalues of M then p {M) = Xr. Since A4" is a non-negative matrix, it h as at least one 
larges t no n-negat i ve eig envalue which coincides with its spectral radius (see iGantmachej 



(12005) or Mevei l 2000h') . When the largest eigenvalue is positive we shall call it, follow 
ing lKimmel and Axelrod j2002l) . the m althusian parameter m = p (M) = of the branch- 
ing process (see also lJagers et a 3 (120071) 1 

The main classification result of multitype branching processes, in the positively regular 
case, s tates that there are only three possible regimes (see lHarrisI (Il963h or lAthreva and NevI 
( 119721 )): sub-critical (m < 1), super-critical {m > 1) and critical (m = 1). Note however, 
that when b = 0, the corresponding branching process is not positively regular. Neverthe- 
less, there is a gen eralization of t he classificati on of multitype branching processes, due 
to Sevastyanov (see lHarrisI (Il963h : ljirinal (Il970l) ). that allows us to include the case b = 
within the same three regimes as above. 



3 Malthusian Parameter of the Phenotypic Model 

In this section we shall employ the perturbation theory for eigenvalues and eigenvectors in 
order to study the average dynamics of the phenotypic model, with b small. In order to main- 
tain the main flow of the text we have omitted most of the computations and esti mates which 
have b ee n placed in the appendix for th e convenience of the reader. See also Iwilkinsoni 
d 19651) or lThompson and McBrid9 d 19741) for more details. 

Metho ds of spectral pe rturbation h ave been used in the study of evolution of macro- 
molecules dSwetina and Schuster., 1982l : lThompson and McBridel , l 19741) . where the the start- 
ing point is a diagonal matrix and the perturbation adds other off-diagonal elements. Our 
approach is similar, but our starting point is an upper triangular matrix and the final matrix 
is tridiagonal. 

In order to implement this program it is necessary to know all the eigenvalues of the 
unperturbed matrix and their associated left and right eigenvectors. The idea is to write 
the mean matrix (|4} as sum of a matrix whose spectral problem is exactly solvable with a 
"perturbation matrix" depending on a small parameter, in such a way that the perturbation 
vanishes as the parameter goes to zero. 

Let M be the mean matrix of (|4} the phenotypic model. Then, since c = 1 — /? — t/we 
may write 

M ==Mo-\-bP, 
where Mq is the mean matrix of the simple phenotypi model: 
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and P is the matrix 



(0 d 
Q[\-d) 2d 
2{\-d) 
Mq= Q 



VO 

/O 
0-10 
1-20 
2-3 
3-4 








Rd 

m-d)j 



(5) 





VO 
































...-(/?- 


-1)0 


... {R- 


1) 0/ 



(6) 



It is easy to see that M is the mean matrix of the phenotypic model when one sets b = 
and c = 1 — (J? in the generating function ([TJ. An interesting feature of this particular case is 
that its spectral problem is completely solvable as functions of the parameters d and R. The 
eigenvalues of the mean matrix Mq are: 



r(l-J) 



:0,...,^. 



In particular, the malthiisian parameter is the largest positive eigenvalue 

= l'^=Rc = R{\-d). 



(7) 



Also, note that the operator norm of P is ||P|| =2(^—1). 

Now we write the eigenvalue X,- of Af as a function of the parameter b, expanded as a 
power series 

A, = A« + A> + A>2 + ..., 

where is the corresponding eigenvalue of Mq. Clearly, — )■ X'j! as Z? — )• 0. The higher 
order perturbation terms A/, {i = 1,2,3,...) are written in terms of all the eigenvalues 
= 0, ...,/?) of Mq and their associated left and right normalized eigenvectors dJ* and uj. 
We are interested in the malthusian parameter m of the matrix M: 

m = irP + m^b + n'P'b^ + . . . , 

where = ( 1 — J)^. A lengthy calculation gives the second order expansion of the malthu- 
sian parameter m: 

\)d 



m={\-d)R + 



{R- 



1 



-Rb 



(l+Jj?/2)(j?- l)^J„,2 , ^,,^3, 
(1-^)3 



(8) 



-Rb^ + 0{b^ 



The parameter space of the model is {{b.d.R) 6 x N | b + d ^ 1}. If we forget the 
discrete parameter R, then this set can be identified with the triangle {(0,0), (0, 1), (1,0)} in 
the (cf,Z))-plane (see ¥IGM- 
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1/2 



2/3 



3/4 



R=4 



Fig. 1 Parameter space (d,b) of the phenotypic model with the critical curves m = 1 for R = 2,3,4. 

Now we can make some quantitative predictions about the behaviour of the model. 

One can localize in the parameter space {b,d) each of the three regimes of the phe- 
notypic model {sub-critical, super-critical, critical) by considering, for each fixed R, the 
critical curve m(b, d,R) = I (see FICfTJ. The region of sub-criticality (for each fixed R) lies 
on the right side of the critical curve and the region of super-criticality (for each fixed R) 
lies on the left side of the critical curve. 

Another interesting feature of the parameter space is related to the points of intersection 
of the critical curve m(b,d,R) = 1 with the boundary curves b->rd = I and b = 0. The points 
of intersection with the boundary curve b = are give by dc{R) = 1 — 1//? and the points 
of intersection with the boundary curve b -\-d = I have non-zero coordinates and can be 
written as (1 —bc,bc) where bc(R) is the ^-coordinate and 1 —bc{R) is the rf-coordinate. For 
instance, it is easy to find the exact value of b,- for R = 2: bc{2) = \ — = 0.2928. This 
calculation gives an important prediction: if b > 0.2929 then the population can not become 
extinct by increasing d towards its maximum value (= 0.7070), with probability going to 1 
as R becomes large. 

More generally, one has the following criterion of no sure extinction. Consider the phe- 
notypic model with parameters {b,d,R), with R large enough. If 7^ is sufficiently small 
and (up to order 0{b^)) 



then, asymptotically almost surely, the population can not become extinct by increasing the 
deleterious probability towards its maximum value. In other words, when the inequality (|9) 
is satisfied the process can have only one regime, namely, the supercritical, which renders 
the model completely trivial, since there is no transition. 



bR{R- \f > 1 



(9) 
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This criterion follows from the approximation of bc(R) that can be obtained from the 
truncated perturbation expansion of the malthusian parameter. Solving the equation m{b, 1 — 
b,R) = I for b in terms of R gives 

Here we should assume that R is sufficiently large since the truncated expansion is a good 
approximation for in only when b is small. In fact, 7? > 10 is reasonable, since dc{lO) = 0.9 
and fof(lO) can be numerically evaluated to give bc{lO} = 0.00486. Using the approximation 
GOll one obtains bc{\0) w 0.00123. 

From now on we shall split the analysis of the phenotypic model according to which it 
is sub-critical, super-critical or critical. 



3.1 The Sub-critical Regime: Lethal Mutagenesis 

The first consequence of our results is a proof, in the context of the phenotypic model (pro- 
vided one assu mes that al l effec ts are of pure ly mutational nat ure), of the conjecture of lethal 
mutagenesis of lSull et~a 1 ( l2007h . Recall that lBuU et aj ilOOlj) assumes that all mutations are 
either neutral or deleterious and write the mutation rate U = Ud + Uc where the component 
Uc comprises the purely neutral mutations and the component comprises the mutations 
with a deleterious fitness effect. Let ^max denote the maximum reproductiv e capacity among 
all particles in the viral population. The extinction criterion proposed by IBuII et all ( |2007|) 
states that a sufficient condition for extinction is 

e"^"^max<l. (11) 

According to lBull et"a i ll200ll2008h . s-'^'i is both the mean fitness level and also the fraction 
of offspring with no non-neutral mutations. In the absence of beneficial mutations the only 
type of non-neutral mutations are the deleterious mutations and hence e^^'' = c = I — d. 
Since the evolution of the mean values depends only on the expected values of the replication 
distribution f,., it follows that Rm^x = R- That is, the extinction criterion (lllb is equivalent, in 
the context of the phenotypic model, to 

{\-d)R<i, (12) 

which is exactly the condition for the model to be sub-critical. In lBull et~a ]( l2008h the authors 
suggest a modification of the extinction threshold dllb that accounts for beneficial effects as 
long as they do not couple the deleterious ones. Assuming that the population size is large 
enough and that a large number of individuals experience the full set of beneficial mutations 
they propose the following threshold 

e-""R^,Al+b)<l. (13) 

Our formula ([Sj for the malthusian parameter allows us to propose a generalization of the 
extinction threshold dllb without any extra assumptions: if Z? ^ is sufficiently small (up to 
order 0{b^)) and 

(l-J)/?(^l + ^^fe) <1 (14) 
then, with probability one, the population become extinct in finite time. 
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Observe that, when b = and R fixed, the model becomes critical for d = ddR) = 
l — l/R.ln other words, for every fixed mean replication capacity R, if the deleterious prob- 
ability satisfies d > dc{R) then lethal mutagenesis occurs almost surely. Thus one may con- 
sider a general critical probability function dc{b,R) and in order to compute this probability 
as a function of b and R one must solve the equation m{b,d,R) = 1 for d in terms of {b,R). 
This can be done by implicit differentiation of the equation m{b,d,R) = 1: 

d,{b,R) = (1 - l/R)+b{R-\f + 0{b^) 

= d,{R)+b{R-lf + 0{b^). 

Formula il5l allows us to obtain a weak generalization of the extinction criterion to the case 
where b is non-zero and d w dc{R). If ^ is sufficiently small and (up to order 0{b^)) 

{\-d)R<\-bR{R-\f (16) 

then, with probability one, the population become extinct in finite time. We summarize the 
discussion about extinction criteria in Table [T] 



Extinction Criterion 



Conditions 



Comments 



(l-^/)R< 1 
(no beneficial mutations) 

(\-d)R{\ + b) 
(uniform beneficial mutations) 

[R- l)d 



[\-d)R 1 + 



(l-rf)2 
(strong version) 



< 1 



(l-fl')R< \-bR(R-\f 
(weak version) 



N ^oo 



The lethal mutagenesis criterion (LMC) 
of BuUetal (2007) with (l-d)= 

Generalization of the LMC proposed in 
Bull et al (2008) (N =population size) 

Branching process generalization 
obtained from perturbative expansion 

Branching process generalization 
solving implict equation m{b,d,R) = 1 



Table 1 Extinction criterion and its generaliations. The role of beneficial mutations h in several forms. 



The main conclusion here is that the existence of lethal mutagenesis depends on "ge- 
netic components" (mutational rates) and other additional deleterious effects (host driven 
pressures intensifications), as well as on strict "ecological components", namely, the maxi- 
mum replication capacity of the particles in the population and on the initial population size. 
As a result the viral population may reach extinction by increasing the number of deleteri- 
ous mutations per replication cycle, by decreasing the value of R in the population or by 
a combination of the two mechanisms. The m utational strategy is the basis of treatments 
using mutagenic drugs (see lCrottv et all (l200lh ) that induce errors in the generation process 
of new viral particles reducing their replication capacity. A straightforward consequence 
of formulas ( 112b . (114b and ( 116b is that a single particle showing the maximum replication 
capacity R is able to rescue a viral population driven to extinction by mutagenic drags. If 
it is assumed that RNA viras populations correspond to a swarm of variants with distinct 
replication capacities, for a therapy to become effective it is important that it will eliminate 
the classes represented by particles with highest replication capacities. In other words, the 
higher the replication capacity of the first particles infecting the organism the larger should 
be the number of deleterious mutations (or effects), that is, the larger should be the drag 
concentration. This can be a clear limitation for treatments based on mutagenic drags. 
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3.2 The Super-critical Case: Relaxation and Equilibrium 

In the super-critical regime, the population grows at a geometric pace indefinitely. Never- 
theless, there are two distinct phases that occur during this growth: a transient phase ("re- 
laxation"or "recovery time") and a dynamical stationary phase. 

3.2.1 Relaxation towards equilibrium 

An important question concerning the adaptation process of a viral population to the host 
environment is the typical time needed to achieve the equilibrium state. As the equilibrium is 
characterized by constant mean replication capacity an obvious criteria to measure the time 
to achieve equili brium would be by the vanishing of the variation of this variable as used 
in other studies jAguirre et all l2009h . Neveitheless, this method is clearly subjected to the 
limitations of numerical accuracy with evident drawbacks if one wants a sharp and universal 
criterion to differentiate populations from the point of view of how fast a population can be 
typically stabilized in a organism. 

Viral populations are commonly submitted to transient regimes. As pointed out earlier 
the infection transmission process represents the passage of a small number of particles from 
one organism to another in such a way that in this process the viral population is submitted 
to a subsequent bottle-neck ejfect during spreading of viruses in the host population. In order 
to approach the problem of relaxation after a bottleneck process in a more sound basis the 
natural quantity to be considered is the characteiistic time derived from the decay of the 
mean auto-correlation function. When the temporal mean correlation function C{n) is of the 
form exp(— na) the decay rate is defined as the parameter a. The natural way to define the 
characteristic time T to achieve equilibrium is by setting 7 = I /a. 

The temporal mean correlation function C{n) may be calculated by considering recur- 
sive applications of the mean matrix M on the initial population: ZqM" Zq. Since we are 
looking for bounds on the mean correlation function, it is enough to co nsider the can onical 
ini tial populations Zq = e^- and use the Perron-Frobenius theorem (see dHarrisl Il963l p. 38) 
or jAthreva and Nevll 19721. p. 185)) to write 

ZlM"Zo = m"{uR{r)vR{r))+0{p'^) , 

where ur is the normalized right eigenvector and vr is normalized left eigenvector (see the 
appendix for more details) corresponding to m and < po < m. Moreover, since we are 

0(p") o(p") 

assuming that m > 1, it follows that < < I and ,,° — )• as « goes to infinity. 

Define the type r of an initial population Zq as the largest replicative class r represented 
in that population. In general, the mean correlation function will depend only on the type 
r of the initial population and will be denoted by Cr{n). Therefore, the mean correlation 
function is typically exponential and given by 

Cr{n) = Kr{n) exp ( — nlog(m)) , 

where K,{n) {r = 0,...,R) is given by 

1 



Mn) = 



UR{r)vR{r} + -J- 



It is difficult, in general, to calculate an explicit expression of UR{r)vR{r), but it is easy to 
see that 

1 , , 1 
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and so the time dependence of Kr is negligible. Indeed, for n sufficiently large the dominating 
factor of C, (h) is the exponential. 

The lower bound for the mean correlation is attained when the type of the initial popu- 
lation is R, that is, 

KRK\ + 0{b) w 1. 
In this case the mean correlation function is 

Cr(«) w exp(-«log(m)) , 

very similar to the case with no beneficial effects. In fact, in this particular case b = and 
one has K, = 1, for every r, so it is enough to consider the case Zq = e^, that is, 

C(«) =exp(-«log(^(l-d))). 

The difference between these two cases lies on the value of the largest eigenvalue (the corre- 
lation decay rate) m. As m{b) > m(0) for the same values of d and R the correlation function 
decays faster if beneficial effects are present. 
In general, we have that (see the appendix) 

KrXO{b-^'^-'''>) for r=l,...,^. 

Hence the mean correlation function may be written as 

'^'■(") ~ ]^ ( - «log(m)) . 

where K,- = K, {d,R). 

The conclusion here is that when the type of the initial population is r there is a "delay 
effect" on the decay of the mean correlation of order relative to the fastest decay 

(lower bound) Cr(w). The slowest possible decay (upper bound) is attained when the type of 
the population is r = 1, which have a delay of order Z?^'*^''. This "delay effect" on the mean 
correlation appears in simulations as "jumps" of magnitude (1 — J) of the mean replicative 
capacity, somewhat reminiscent of the "punctuated equilibrium scenario". 

Observe that the closer is the parameter d to its critical value dc, the longer is the time 
needed to achieve equilibrium. The clearest way to characterize the time behavior of the 
viral population at or around the critical point is through the typical time T to approach 
equilibrium derived from the decay of correlations described above. The expression T = 
1 / log (m) shows that at the critical point the equilibrium state is never reached, i.e. , the decay 
to equilibrium is at least non-exponential. A scaling exponent characterizing the behavior of 
T in the neighborhood of the critical point dc can be easily obtained. The expansion around 
dc gives 

m'[dc) 

In the particular case where = we have 

TK{l-dc)\d-dc\'\ 

and hence l/m'{dc) = {l-dc) + 0{b). 

One possible application of this result relates to the very initial phase of the infection 
process. If one considers that during this phase the host immune system has not been yet 
stimulated against the vims, one can assume that the deleterious effects would be solely 
represented by the viral intrinsic mutation rates. Therefore, the largest the value of R, i.e., 
the largest the replication capacity of the initial viral particle, the fastest the decay of the 
progeny auto-correlation. Intuitively the parameter R defines the degree of virulence of the 
infection during the early stage of the infective process. 
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3.2.2 The Dynamical Stationary State 

According to the "Malthusian Law of Growth" it is expected that a super-critical branching 
process will grow indefinitely at a geometric rate proportional to m", that is, Z„ w m" W„, 
where W„ is a random vector. That is, the normalized random vector W„ = Z„/m" posses 
a finite limit when n —5- °°. More precisely, there exists a random variable W ^0 such that, 
with probability 1, 

\\mW„ = W. 

Moreover, W = Wu where u is the normalized right eigenvector corresponding to the 
malthusian parameter m and 

E{W\Zo) = {v,Zo), 

where v is the left eigenvector corresponding to the malthusian parameter m. The meaning 
of this result is that the total size of the population divided by ?n" , converges to a random 
vector, but the rel ative frequencies of the classes of particles approach fixed limits. In fact, 
( iKurtz et all 1 1994 ) formalized this statement as a limit theorem (no averages here!): 

Zyi 

lim r = II (almost surely). (17) 

|Z„ 

This result is useful in computational simulations of the model, since one may find an ap- 
proximation to the "deterministic part" of W by sampling the population and computing the 
relative frequencies of replicative classes. 

Recall that the normalized right eigenvector = (i(^(0), . . . , mr(^)) corresponding to 
the malthusian parameter m is positive and is normalized so that Lr"R('') = therefore, 
it may be seen as a probability distribution on the set of classes {0, . . . ,^}. It is called the 
mean fitness distribution of the replicative classes of the viral population. 

Using the same perturbative techniques as before, it is possible to compute the perturba- 
tion expansion of the right eigenvector corresponding to the malthusian parameter. We are 
interested in the right eigenvector iis associated to the malthusian parameter m. The first 
order expansion of the normalized right eigenvector ur associated to m may be written as 

/t=0 

where the terms jS^; are functions of d and R. The explicit expressions of the terms jSj- are 
somewhat cumbersome, however we still can use the abstract formula to gain some insight 
about ur. 

The normalized right eigenvector of Mq is 

u\{r) = binom(r;/?, \-d) = {1){\- d)' d'^-' . 

From these formulas one immediately obtains, for b = 0: 

- The mean replication capacity E(tx^) = R{\ — d). 

— The phenotypic diversity Var(u^) = Rd[ \ — d). 

In general, the mean replication capacity of a viral population is given by 



R 

k=0 



(19) 
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Indeed, the malthusian parameter may be calculated as 

KZ„+i)| \M{Z„)\ 
m = lim — — — = lim ——, — -— . 

\{Z„)\ n^oo |(Z„)| 

Now it is easy to see that \M{Z„) \ = Y,kkZ^ and hence by the Perron-Frobenius theorem, 
equation i\% follows immediately. In particular, when is small enough: 

E{uR)=R{i-d) + ^^f^Rb + 0(b^). 
L — a 

Using ( 118b one may write the phenotypic diversity of a viral population as: 

Var{u^)=/?^/(I-J)+|^2^ftCov(u^,u°)j b + 0{b^). (20) 

It is possible to estimate the magnitude of the first order perturbation term in eq. (I20t : 

Var(..) « /..(I - (l + 8.^ ^^^^i±(i^) . 

It is well accepted that the phenotypic diversity is an important c haracteristic of the 
viral population i ntuitively related to the idea of population robustness dForsteretallbOOd 
lElena et aiLl2007h . In fact, a homogeneous population would be less flexible from the point 
of view of adaptation. The variance associated with the stationary state can be seen as a 
measure of phenotypic diversity. It shows that, in the case where b = 0, the maximum value 
of the phenotypic diversity R/4is reached if d = 1/2 for any value of R. For b =^Oit can be 
shown numerically, that there is always a value Jmax of d for which the population attains 
maximal phenotypic diversity. This value d^^^^ is a decreasing function of b. For w this 
value is Jniax ^ 1 /2. 

It is interesting t o note that experimental me asurements of c = 1 — t/ are close to 0.5 
dSaniuan et alll2004l : iDomingo-Calap et alll200^ . suggesting that viral populations follow 
a principle of maximal phenotypic diversity. 

As far as the experimental situation of the phenotypic diversity is con cerned, one of th e 
first attempts to experimentally measure the phenotypic diversity was in iDelbriickl IT945I) . 
where the total "burst size" of a progeny from phage-infected bacterial cells was estimated. 
More recently, measurements of th e phenotypic diversity in vitro generat ed by single viral 
particles infecting individual cells jZhu et all l2009l : Irimm and Yinl l2012l) revealed a broad 
distribution of virus yields (50 to 8000 progeny virus particles). One may regard these re- 
sults as indication that the replicative capacity of a virus from a particular replicative class 
is characterized by a fitness distribution obtained as the distribution of progeny produced 
by representative particles from that same replicative class. Although these authors did not 
went further and investigated if different particles form a viral population have different fit- 
ness distributions their results suggest that a fitness distribution over a viral population may 
resemble the mean fitness distribution of the replicative classes obtained from the pheno- 
typic model. In this direction, further carefully designed experiments aiming at determine 
the mean fitness distribution of a viral population would be welcome. 
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3.3 Criticality and Eigen's Error Threshold 

The relation betx yeen multity pe branching processes and Eigen's theory of evolution of 
macromolecules (lEigenLll97lh has been established bv lDemetrius et all ( Il985l) . In this work 
it is shown that one may canonically associate a certain difference equation to a discrete time 
multitype branching process through its mean matrix. This deterministic dynamical system 
is called Eigen 's selection equation, due to its similarity to the phenomenological equation 
describing the kinetics of self-reproducing molecules in a dialysis reactor. Moreover, the se- 
lection equation is essentially equivalent to a linear difference equation x,, = M"xo, where 
M is the mean matrix of the branching process. Note that this equation is formally identi- 
cal to the equation for the evolution of the mean values of branching process (see eq. ((3)). 
Hence the selection equation may be seen as a mean field limit equation associated branch- 
ing process. In particular, when the process does not become extinct, part of its asymptotic 
behaviour, namely, the relative frequencies the classes, can be recovered by the selection 
equation. Moreover, it is shown that a super-critical branching process displays "freezing 
in" of initial fluctuations, that is, the coefficient of variation of the process vanishes asymp- 
totically almost sure. In this sense, if the population is infinite and one waits long enough 
before starting observation, the deterministic model is fairly reliable because fluctuations 
are small. Nonetheless, when considering finite populations, i.e., finite size samples of each 
generation, the deterministic approximation is only part of the story, since the fluctuations 
contain the "out-of-equilibrium" characteristics of the system. 

The use of branching processes links th e concept of error thre shold with that of criti- 
cality. If we switch to the genotypic view of IPemetrius et all ll 19851^ and suppose that there 
are no back-mutations (b = 0) then we can consider polymers with chain length of v nu- 
cleotides and that there is a fixed probability p for copying a single nucleotide correctly, 
that is, c = . The formula for the malthusian parameter (m = cR = p^R) together with the 
criticality condition {m= 1) gives the stochastic error threshold of jPemetrius et alL[l985L 
p. 254, eq. (47)): 

InR 

(21) 



-Inp 



Since the lethal mutagenesis criterion, in the context of branching processes, is exactly the 
criticality condition, the stochastic error threshold refers to the probability of extinction. This 
can be compared with the determinis tic error threshold (see .Eigea (,1971,) ; ,Eigen and Schusteil 
l ll979l) : ISwetina and Schuste^ ( Il982h ): 

v, = ^, (22) 
— Inp 

where the parameter a is the superiority of the master sequence (see lDemetrius et a lliiil). 
The deterministic error threshold is based on the condition that the error-free productivity 
of the master molecular species becomes equal to the mean total productivity of all other 
species. The condition to replicate with a fidelity above the error threshold will always be 
valid for the master species provided the mutation terms of all other species are sufficiently 
small. Beyond that point, the master species is no longer maintained: the population has 
reached the error catastrophe (Bull et al, 2005). Moreover, the demand to operate above the 
stochastic threshold is always a stronger condition than the corresponding requirement of 
the deterministic equation. 
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It is interesting to note that at the critical point the dynamics of the branching process 
relax to equilibrium according to the typical algebraic decay in time. This property raises the 
question (unapproachable in the present context) if the classes in the critical population are 
so strongly correlated in such a way that the whole population behaves as if it is composed 
be only one replicative class. This behaviour reminds one of the basic p roperties of the error 
threshold in Eigen's theory of macromolecular evolution (lEigenLll97lh. 

Prelimina ry results concerning the dynamics of fluctuations jAntoneli et a 1 ( l2013bl) : 
ICastrol ( I2OIII) ) indeed show that the time variation of the numbers of particles in each sep- 
arated class follows a pattern such that the variation observed in one class is rigorously the 
same observed in all the others. This indicate a high level of correlation between the classes 
in complete analogy with critical phenomena of many physical systems. It is worth to note 
that, in fact, there is a correspondence between Eigen's m odel and the equilib rium statistical 
mechanics of a certain inhomogeneous Ising system (see lLeuthaussej (ll987h V 



4 Outlook 

Modeling the dynamics of viral populations by means of multivariate branching processes 
does not take into account many molecular/microscopic details of the replicative process, 
nevertheless, it is remarkable that the provided description at the population/macroscopic 
scale is very appealing and well fitted to various aspects of the observed phenomenology 
of viral systems. This leads to the conclusion that we are probably far from exhaust the 
analytical and predictive power of branching processes as a mathematical tool for the study 
of viral populations. In fact, many aspects of the problem like (i) the role of finite size 
effects, (ii) the role of different infected cell types for the problem of evolution of viral 
populations in multi-cellular hosts, (iii) the relation between quiescent infected cells with 
branching process with singular components and (iv) the descriptive limits of discrete versus 
continuous time branching processes, are still to be considered by future investigations. 
Also it is worth to mention that most of the quantitative results and all of the qualitative 
observatio ns described he re can be easily reproduced by computer s i mulat ion of the model 
(see Castro ICastrol (l201lh and manuscript in preparation ICastro et~al ( l201lh '). 

Finally, it is important to emphasize that the branching process analysis carried out here 
revolves around the malthusian parameter, which is based on certain critical assumptions re- 
garding environmental constraints, for instance, that the environment (host) does not change 
dramatically. The significance of this macroscopic m easure in stu d ies of evolutionary pro- 
cesses has been challenged in a seri es of articles (see lDiet3 ( l2005l) ; lKowald and Demetrius! 
(l2005h ; IZiehe and DemetriusI ( |2005|) ). When the environment is strongly perturbed, another 
parameters become important, as the "evolutionary entropy", which is a measure of the vari- 
ability i n the age at which individual s produce offspring and die (see lKowald and DemetriusI 
1 I2OO5I) : IZiehe and DemetriusI ( l2005h ). These parameters will play a role in studies of the 
extinction process. 
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Appendix: Spectral Perturbation Theory 

In order to carry out the perturbative calculations we need to solve the spectral problem for the matrix 



Mq-. 



(0 d 
Q (l-d) 2d 
2(1 -rf) 




\ 




Rd 
R(\-d)J 



\0 
The eigenvalues of the mean matrix Mq are: 

X'^ = rc = r(\ - d) r = 0,...,R. 
In particular, the malthusian parameter is the largest positive eigenvalue 

m° = X'^=Rc = R(\-d). 
The normalized left and right eigenvectors, uj? and tt[l, associated to the eigenvalue A^", satisfy: 



Wilting and it J? in components as 



^(v?(0),v«(l)„ 
^K(0),«°(1), 



one has that 





-d{r+\) 
\ (1 

{\-df 



for 
for 



k = 0,...,r-\, 
k = r+\ 



for k = r,r + 2,... ,R. 



and 



Qil-dfd''-'^ for k = Q,...,r, 
for k = r+l,....R. 

Now we write the eigenvalue A,, of M as a function of the parameter b, expanded as a power series 

X,- = X° + X}h + X^h^ + ... , 

where A,!* is the corresponding eigenvalue of Mq. Clearly, Xi-{b) — > A'' as ii — > 0. The higher order perturbation 
coefficients A/, (; = 1, 2,3, . . .) are written in terms of all the eigenvalues A," (s = 0, . . . ,R) of Mq and their 
associated left and right normalized eigenvectors i;^ and and the perturbation matrix 



/O 
0-10 
1-20 
2-3 
3 -4 




Vo 



-(fi-1) 
(R-1) O/ 



Note that the operator norm of P is \\P\\ = 2{R — I) and so the magnitude of the perturbation is 2b{R — 1) 
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The general expressions for the first and second order coefficients of the perturbation expansion of an 
eigenvalue are 



We want to use these formulas to compute a perturbation approximation (for b around 0) of the malthusian 
parameter m(b) of (M). Recall that, m(0) = m" = S(l — d) and the corresponding left and right eigenvectors 
are given by 

^» = 1/(1-<(0,...,0,1), 
u° = {uk(0),...,uk{R)), 

with u/i(A) = binom{k;R, l—d). The first order coefficient of the malthusian parameter is 

For the second order coefficient we also need 

= l/(l-rf)'^-'(0,...,0,l,-Rrf/(l-£/)). 

The second order coefficient is 



„2 _ 



5 5 
/=0 '^R ^ '^i 



R{R-\f{\ + ^d)d 

Analogous perturbative formulas exist for the perturbation expansion of the left and light eigenvector 
con'esponding to an eigenvalue A^. The left and right eigenvector i>, and it,, of the matrix M associated with 
the eigenvalue A, can be written as a function of the parameter b, expanded as a vector-valued power series: 

= + £ bIu^ + 1? £ B;. ^^0 ^ ^ ^ ^^ 



where v]! and ttj: are, respectively, the left and right eigenvectors of the matrix Mo associated to the eigen- 
value Xj}. The perturbation terms A[ , and with i = 1,2,3,..., can be written in terms of the eigenvalues 
of Mo and their associated left and light eigenvectors and the perturbation matrix P. Observe that when 
— > we get the normalized left and right eigenvectors and tt^ associated to the dominant eigenvalue 
m" = (1 — d)R of the mean matrix Mo- 

We are interested in the normalized eigenvectors vn and ur associated to the malthusian parameter 
m = {l-d)R 

VR = vl + bj^a,vl + 0(J?), 

k=0 

with at = aJ; J, (0 s; A: s; R - 1) and = - ifl^' aj.. 



u„ = ul + bf^p,u° + 0{b^), 

k=0 
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with lit = Bj; J., (0 s: <: s; R - 1 ) and = - ifl^ ft 
The first order coefficients are given by 



AO -A? 



and the second order ternis are given by 



fo (A«"-AP)(A^"-A,") 



(AO -AO) (AO -AO) 



fori = 0,...,R- 1. 

One may use the fonnula of up to second order 

(i?-l) 



(i-ciy 



and the first order expansion = tt^ + 0{h) in order to estimate the product of coefficients iiR(r)VR(r) (up 
to their leading order). For instance, it is easy to find that 

UR{R)vR(R) = \ + 0{h) = 0{\), 
UR{R-l)vR(R-\) = 0[h).. 
ur{R-2)vr(R-2) = 0[1?). 

In general, using a higher order expansion formula for vr, one concludes that 

itR(R-r)vR{R-r) = 0(¥). 

This follows form the fact that the coefficient of order ¥ in the perturbation expansion of vr is a linear 
combination of the left eigenvectors , . . . jijO^^, all of which have zeros in the first r — 1 components. 
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